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ABSTRACT 

We present the results of two-dimensional and three-dimensional magnetohydrody- 
namical numerical simulations of relativistic magnetic reconnection, with particular 
emphasis on the dynamics of the plasma in a Petschek-type configuration with high 
Lundquist numbers, S ~ 10 5 — 10 s . The numerical scheme adopted, allowing for 
unprecedented accuracy for this type of calculations, is based on high order finite vol- 
ume a nd discontinuous Galcrkin methods as recently proposed bv lDumbser &: Zanottil 
(2008). The possibility of producing high Lorentz factors is discussed, showing that 
Lorentz factors close to ~ 4 can be produced for a plasma parameter a m = 20. More- 
over, we find that the Sweet-Parker layers are unstable, generating secondary magnetic 
islands, but only for S > S c ~ 10 s , much larger than what is reported in the New- 
tonian regime. Finally, the effects of a mildly anisotropic Ohm law are considered in 
a configuration with a guide magnetic field. Such effects produce only slightly faster 
reconnection rates and Lorentz factors of about 1% larger with respect to the perfectly 
isotropic Ohm law. 
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1 INTRODUCTION 

Relativistic magnetic reconnection is recognized to be a 
key physical process in high-energy astrophysics, being able 
to convert magnetic energy into heat and plasma kinetic 
energy over short timescales. Relevant examples include: 

(i) the magnetospheres of pulsars near the Y-point, where 
the outer most magnetic field lines interse ct the equato- 
rial plane (Uzdenskv 20031; iGruzinovl 12005); (ii) the dissi- 
pation of alternating fields at t he termination shock of a 
relativistic striped pulsar wind (|Petri fe Lvubarskvl l2007l h 

(ii) soft gamma-ray repeaters, where giant magnetic flares 
could be the explanation of the observe d strongly mag- 
netiz ed and relativistic ejection events (jLvutikovl |2003| . 
2006); (iv) gamma-ray burst jets, where particle acceler- 
ation by magnetic reconnectio n in electron-positron plas- 
mas is supposed to take place (ID rcnkhahn fc Sprufq 2002; 
Barkov fc Komissarovl l20ld ; iMcKinnev fc Uzdenskv! hold ; 

(v) and accretion disc coronae of ac- 
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tive galactic nuclei, where violent releases of energy may be 
generated by Petschek magnetic reconnection of strong mag- 
netic loops emerging from the disc via buoyancy instabil- 
ity (|di Matted 1 19981 ; ISchopper et~ailll998l ; IJaroschek et all 



* E-mail: zanotti@aei.mpg.de 

1 A similar mechanism has been proposed by iTanuma et al.l 
to explain the X-ray emission in the Galactic halo. 



12004 ). Very recently, moreover. iNalewaiko et all |201ll ) have 
described a model of mini jets in blazars to explain their 
ultra-fast variability and were able to fit data of PKS2155- 
304 assuming that the dynamics is governed by relativistic 
magnetic reconnection with a weak "guide magnetic field", 
i.e. with a magnetic field aligned to the current. 

In the recent past, both theoretical studies and numer- 
ical investigations have greatly improved our understanding 
of relativistic magnetic recon nection. Although within the 
incompressibility assumption, iLvutikov fc Uzdenskv! (|2003l ) 
first found that in the Sweet-Parker reconnection three dif- 
ferent regimes can be produced, which depend on the ratio 
of the magnetization pa rameter a m to the Lundquist num- 
ber S. iLvubarskvl ([2005), on the other hand, explicitly ad- 
dressed the question about whether the reconnection rate 
can be significantly enhanced or not in the transition from 
Newtonian to relativistic magnetic reconnection. He found 
that, unless the reconnecting fields are not strictly antipar- 
allel, the relativistic Petschek reconnection should not be 
considered as a mechanism for the direct conversion of the 
magnetic energy into the plasma energy and the reconnec- 
tion rate would be at most 0.1 t he speed of light, contrar y 
to what o riginally sugges t ed by iBlackman fc Field! |l994l ). 
Moreover, iTolstvkh et al.1 |2007l ) extended the Petschek re- 
connection model to incorporate relativistic effects of im- 
pulsive reconnection and, for current layers embedded into 



2 0. Zanotti, M. Dumbser 



strong magnetic fields, they claimed that the plasma can be 
accelerated to high Lorentz factors. 

Because our understanding of relativistic magnetic re- 
connection is still incomplete, in the last few years there has 
been a growing expectation towards numerical simulations 
as a promising tool for clarifying the rich underlying physics, 
particularly in the non-linear regi me. After the pioneering 
resisti ve relativistic simulations by IWatanabe fc Yokovamal 
(2006), who considered the Petschek type reconnection in 
the relativistic regi me with a res i stive m agnetohydrody- 
namic (MHD) code. lZenitani et all (|2009bh investigated the 
relativistic reconnection in an electron-positron plasma by 
a two-fluid MHD code and showed that the reconnection 
rate is higher and high er in magnetically dominated regimes. 
IZenitani et al.l (|2009al ). on the other hand, clarified the role 
of a guide magnetic field, which is essentially that of mak- 
ing the output energy flux Poynting dominated rath e r than 
enthalpy dominated. Very recently, IZenitani et al.l (|2010l ) 
showed that, when the resistivity is current dependent, plas- 
moids are repeatedly formed in the current sheet. In spite 
of these progresses, two major numerical limitations still 
prevent the application of numerical calculations to real- 
istic physical and astrophysical scenarios in the relativistic 
regime. The first limitation manifests when treating plas- 
mas with high magnetizations parameters a m , while the sec- 
ond limitation manifests when high Lundquist numbers S 
are encountered, being related to the stiffness of the equa- 
tions in this regime. Physically, high Lundquist number plas- 
mas are very interesting as they are supposed to become 
unstable and break up into secondary magnetic islands. 
In the Newtonian framework, for example, ISamtanev et al.l 
(|20ogT ) showed that Sweet-Parker current sheets are unsta- 
ble to super- Alfvenically fast formation of plasmoid chains 
for 3 > S c ~ 10 4 , and it is not clear the extent to which 
this instability persists in a relativistic framework. 

In this paper, by adopting the inn o vative numerical 
method presented in lDumbser fc Zanottil (|2009h . we recon- 
sider some of the scenarios discussed above within the rel- 
ativistic magnetohydrodynamic single fluid approximation, 
focusing, in particular, on the dynamics of the plasma in 
a regime characterized by high Lundquist numbers, S ~ 
10 5 - 10 8 and mild magnetization parameters a m ~ 1 — 20. 

We also investigate numerically, for the first time, the 
effect of an anisotropic Ohm law, i.e. when the conductivity 
of the plasma is a tensor that depends on the orientation of 
the magnetic field. 

The plan of the paper is the following. In Section [2] we 
report the relativistic resistive equations and the basic phys- 
ical assumptions, while Section [3] is devoted to a succinct 
presentation of the numerical method and to the validation 
of the code in three space dimensions. Section [3] contains 
the results of our analysis, and Section[5]the conclusions. We 
have considered only flat spacetimes in pseudo-Cartesian co- 
ordinates, namely the metric 77 M „ = diag(— 1, 1, 1, 1), where 
from now onwards we agree to use Greek letters p,, v, A, . . . 
(running from to 3) for indices of four-dimensional space- 
time tensors, while using Latin letters i,j,k,... (running 
from 1 to 3) for indices of three-dimensional spatial tensors. 
We set the speed of light c = 1 and make use of the Lorentz- 
Heaviside notation for the electromagnetic quantities, such 
that all V^Lt factors disappear. Finally, we use Einstein sum- 
mation convention over repeated indices. 



2 MATHEMATICAL FORMULATION 

The total energy-momentum tensor of the plasma that we 
consider is made up by two contributions, = + T^ I/ . 
The first one is due to matter 

7T = hpuW+pq*", (1) 

where u M is the four velocity of the fluid, while h, p and p are 
the specific enthalpy, the rest mass density and the pressure 
as measured in the co-moving frame of the fluid. The second 
contribution comes from the electromagnetic field 

Tf = F» x F vX - ^(F XK F XK )r,^ , (2) 

where F M ", and its dual F*' 1 ", is the electromagnetic tensor 
given by 

F^ = n^-BV + e^BA (3) 
F*^ = n"B"-BV- t '"' A ' t £ A n K . (4) 

E v and B" are the electric and magnetic field as measured 
by the observer defined by the four- velocity vector ra M , while 
e y,v\K _ ^ V \ K ^ j s £ ne completely antisymmetric spacetime 
Levi-Civita tensor, with the convention that e 0123 = 1. If W e 
now set n M to define the inertial laboratory observer, namely 
n M = (1,0,0,0), normalized such that n M n M = —1, then the 
four vectors of the electric and of the magnetic field are 
purely spatial, i.e. E° = B° = 0, E i = E t , B l = B t in this 
frame. On the other hand, the fluid four velocity u 11 and the 
standard three velocity in the laboratory frame are related 
as v = v l — u 1 /r, where r = (1 — v 2 )' 1 ^ 2 is the Lorentz 
factor of the fluid with respect to the laboratory frame. 

In Cartesian coordinates, using the abbreviations ft = 
and ft = , the full system of Euler and Maxwell equa- 



tions are 

ftD + ft(D?/) = 0, (5) 

dtSj + diZ] = 0, (6) 

ftr + d,S l = 0, (7) 

d t E i -e iik d j B k + d i V = -J i , (8) 

ftB* + e iik djE h + &<S> = 0, (9) 

ft* + d % E l =p c - ac*, (10) 

ft$ + d t B l = (11) 

dtp c + d i r = 0, (12) 

where the conservative variables of the fluid are 

D = P T, (13) 

S' = wrV + e ijk EjB k , (14) 

t = uT 2 -p+\{E 2 + B 2 ) , (15) 



expressing, respectively, the relativistic mass density, the 
momentum density and the total energy density. The spatial 
tensor Z\ in (J6j , representing the momentum flux density, 
is 

Z) = u;rV Vj - E i Ej - B i B 3 + 

[p+l 2 {E 2 + B 2 )}5), (16) 

where uj — hp is the enthalpy of the fluid while 8] is the 
Kronecker delta. An equation of state is needed to close the 
system, and we have chosen that of an ideal gas, namely 

p=(j-l)pe = ji(u- p), (17) 
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where 7 is the adiabatic index, 71 = (7 — l)/"f, e is the 
specific internal energy. 

While writing equations (|10|l and we have adopted 
the so called dive rgence- cleaning approach presented in 
iDedner et all {2002), which amounts to introducing two ad- 
ditional scalar fields ^ and <I> that propagate away the de- 
viations of the divergences of the electric and of the mag- 
netic fields from the values prescribed by Maxwell's equa- 
ti ons. Additional det ails a bout this approa c h can be found 
in iKomissarovl 1120071 ) and iPalenzuela et af] l|2009h . 

In its most general form, the relativistic forrnula- 
tion of Ohm law is a non- linear propagation equation 
(|Kandus fc Tsagasl [20081, which makes it manifest the con- 
nection between the Ohm law and the equations of mo- 
tion. However, in order to keep the equations numerically 
tractable and because of the poor knowledge of the con- 
ductivity in realistic conditions, simpler forms of the Ohm 
law are usually considered, in which the currents are al- 
gebraically related to the electromagnetic field. Following 
iBekenstein fc Oronl (| 19781 ) we allow for an anisotropic Ohm 
law and therefore we write the four-current vector as 

J M = q u" + ff^c, (18) 

where qo is the charge density in the co-moving frame of the 
plasma, is the conductivity tensor and e„ is the elec- 
tric field in the co-moving frame. Within the collision-time 
approximation the most general form of the conductivity 
tensor is given by 

a"" = a(g^ + fb^V + ^ XK u x b K ) (19) 

where u M is the four velocity of the plasma and 6 M is the mag- 
netic field in the co-moving frame. Th e parameter £ is re- 
lated to the micro-physics of the plasma (|Bekenstein fc Oronl 
1 19781 ) via £ = er/m e , where e and rn e are electron's charge 
and mass, while r is the collision time. As a first applica- 
tion to the case of an anisotropic Ohm law, in this paper we 
consider the case <r MI/ = a{g^ v + £ 2 & M 6"), namely we drop 
the third term in ()19|) . We note, therefore, that anisotropic 
effects in the Ohm law are more and more important as the 
intensity of the magnetic field is increased. The four-current 
vector can also be decomposed in components parallel and 
perpendicular to the observer as 7 M = p c n M + J M , where p c is 
the charge density in the laboratory frame. Hence, it is easy 
to derive from (|18p (see Appendix [S] for the calculations) 
the following expression for the spatial current vector to be 
used in Eq.(JH]), namely 

J = p c v + Ta[E + v x B — (E ■ v)v\ + 

raf(E-B){B~vxE~{B-v)v\. (20) 

Two relevant comments are worth doing about the Ohm law 
(|20|) . The first and most obvious one is that the isotropic 
regime is recovered when £ = 0, in which case E q. (|20[) re- 
duces to the usual expression (jKomissarovl 2007). The sec- 
ond comment is that, in the comoving frame, the Ohm law 
(I20[) becomes 

J = aE + a£, 2 (E-B)B (21) 

which clarifies how the current term proportional to £ 2 is 
present only for configurations for which E ■ B 7^ and it is 
responsible for an extra current term in the direction parallel 
to the magnetic field. 
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Figure 1. Numerical grid of the 3D simulation in the current 
sheet test. The B v component of the magnetic field is reported. 
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Figure 2. Comparison of the numerical solution with the analytic 
one in the 3D current sheet test. The plot shows a one-dimensional 
cut of the magnetic field B y at time t = 10. 

3 NUMERICAL METHOD 
3.1 Brief description 

A well known and challenging feature of the system of equa- 
tions ([5)l- (|12p is that it has source terms in the three equa- 
tions (|8} for the evolution of the electric field that be- 
come stiff in the limit of high conductivity. This pathol- 
ogy has been handled in the recent past by r esorting 
to Strang-splitting techniques JKomissarovl 12(3071) or to 
impli cit-explicit Runge Kutta methods (|Palenzuela et"all 
120091 ). In the present pa per, on the other hand, we apply 
the strategy outlined by iDumbser fc Zanottil |2009), who 
used the so called high order PnPm methods, which are 
a unification of high order finite volume and discontinuous 
Galerkin finite element schemes in a more general fr ame- 
work, and comb i ning r esults from lDumbser et alj l|2008l ) and 
IDumbser et al.l l|2008l ). It is worth stressing that, because 
of the high accuracy they allow to achieve, Galerkin meth- 
ods have been recently considered as a valuable approach 
even in fully relat ivistic calcu l ations . Pro mising results have 
been o btained by IZumbuschl (|2009l ) and iRadice fc Rezzolial 
l|201ll ). 

In our specific implementation, the numerical solution 



4 0. Zanotti, M. Dumbser 



Table 1. Main properties of the initial models. From left to right the columns report the name of the model, the magnetization parameter 
<Tm, the Lundquist number S, the background resistivity r)^, the value of the vertical magnetic field B z , the number of cells used and the 
minimum cell size -L m ; n (used for time-step calculation). All of the models have adiabatic index 7 = 4/3, po = lj PO = 1 and 7]io = 1.0. 
A third order P0P2 finite volume scheme has been adopted. 



Model 


0" m 


5 


Vb 




cells 


-^min 




2D-ml-S1.60e5 


1.0 


1.60 x 10 5 


10" 3 


0.0 


844194 


4.90 x 10" 


3 


2D-m5-S2 . 45e5 


5.0 


2.45 x 10 5 


10" 3 


0.0 


844194 


4.90 x 10" 


-3 


2D-mlO-S2.68e5 


10.0 


2.68 x 10 5 


10" 3 


0.0 


844194 


4.90 x 10" 


-3 


2D-ml5-S2.77e5 


15.0 


2.77 x 10 5 


10~ 3 


0.0 


844194 


4.90 x 10" 


-3 


2D-m20-S2.82e5 


20.0 


2.82 x 10 5 


10" 3 


0.0 


844194 


4.90 x 10" 


-3 


2D-mlO-S2.68e6 


10.0 


2.68 x 10 6 


10" 4 


0.0 


844194 


4.90 x 10" 




2D-mlO-S2 . 68e7 


10.0 


2.68 x 10 7 


10" 5 


0.0 


844194 


4.90 x 10" 


-3 


2D-mlO-S2 . 68e8 


10.0 


2.68 x 10 8 


10" 6 


0.0 


844194 


4.90 x 10" 


-3 


3D-ml . 25-BzO 


1.25 


1.73 x 10 4 


0.005 


0.0 


6105345 


2.49 x 10" 


> 


3D-ml . 25-BzO . 5 


1.25 


1.73 x 10 4 


0.005 


0.5 


6105345 


2.49 x 10" 


-2 



of the vector of conserved quantities is represented at the be- 
ginning of each time-step by polynomials of degree N. How- 
ever, the time evolution of these data and the computation 
of the corresponding numerical fluxes are performed with 
a different set of piecewise polynomials of degree M ^ N, 
which are reconstructed starting from the underlying N de- 
gree polynomials. The part of the algorithm performing the 
time evolution of the reconstructed polynomials uses a lo- 
cal space-time discontinuous Galerkin (DG) finite element 
scheme which provides a local predictor for constructing the 
solution of the par tial differentia l equ ations in the small, 
as it was called by lHarten et all (|l987f l. Our local space- 
time DG predictor allows us to discretize the stiff source 
terms arising in the resistive relativistic MHD equations nat- 
urally and without using any splitting techniques. The local 
space-time predictors are obtained after solving, for each 
individual cell element, a system of non-linear equations, 
which are therefore not globally coupled as in the global 
an d fully implicit space-time Ga l erkin approach introduced 
by van der Vegt fc van der Venl l|2002r ) and by iKlaii et all 
(2006). Once computed, the local space-time predictors are 
then inserted into a global corrector, which is fully explicit 
and provides the coupling between neighboring cells. The re- 
sulting Galerkin scheme can allow for an arbitrary (at least 
in principle) order of convergence, through a one-step and 
quadrature free time update (no need for Runge Kutta time 
stepping) formula. Within this new approach, traditional 
finite volume schemes with N — and usual discontinuous 
Galerkin methods with N = M are included as special cases, 
while the new class of methods with N ^ 0, M > N are in 
general computationally more efficient. In most of the nu- 
merical simulations reported in this paper we have adopted 
the schemes P0P2, P1P2 and P1P3. 



3.2 Validation of the code in three dimensions 



In lDumbser fc Zanottil ((2009j) we presented a wide class of 
numerical tests both in one and in two spatial dimensions, 
showing the ability of the scheme in dealing with the stiff 
terms inherent in the resistive relativistic equations, while 
retaining its high order properties. Here we complete the full 



validation of the code by considering the numerical evolution 
of a self similar current sheet in t hree spatial dimen sions. 
This configuration, first proposed bv lKomissaro vl (|2007l ), has 
the following exact analytical solution for the {/-component 
of the magnetic field: 



By fat) = Bo erf ( ^J~% 



(22) 



where erf is the error function. The initial time for this 
test case is chosen to be t = 1 and the initial condition 
at t = 1 is given by p = 1, p = 50, E — v = and 



B 



(0,B y (x, l),0y . We choose 7 



and Bo = 1. The 



conductivity is set to a = 100, which means a moderate 
resistivity. We have solved the problem with the fourth or- 
der P1P3 scheme on a very coarse unstructured mesh com- 
posed by 3209 tetrahedrons. The grid extension is given by 
[-1.5, 1.5] x [-0.5, 0.5] x [-0.25, 0.25] and it is shown in Fig.|T] 
with a color rendering of the component B y of the magnetic 
field. Fig. [2l on the other hand, shows the perfect matching 
of the numerical solution against the analytic one at time 
t = 10 and computed along a representative one-dimensional 
cut of the numerical domain. We have also verified that the 
method provides the expected order of accuracy when the 
order of the polynomials in the PnPm scheme is changed. 
This concludes the validation of the code in the full three- 
dimensional case. 



4 MAGNETIC RECONNECTION 

4.1 Initial model and boundary conditions 

The initial model that we have considered in our analy- 
sis of relativistic magneti c reconnection is b u ilt on Har- 
ris model, as reported by iKirk fc Skiaeraasenl ((2003) , and 
it reproduces a current sheet configuration in the x — y 
pla ne. Very similar configu r ations have bee n studied also 
by IWatanabe fc Yokovamal |2006) and by IZenitani et alj 
l|2009al lbl). Gas pressure and density are given by p = 
po + CT m po[po cosh 2 (2s)] _1 , p = po + a m p [po cosh 2 (2a;)]" 1 , 
where po and po are the constant values outside the cur- 
rent sheet, whose thickness is 8 = 1. The magnetic field 
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Figure 3. Color map of the rest mass density for two models with the same background electrical resistivity r]b = 10 3 but different 
values of the magnetization: u m = 1.0 (left panel) and <r m = 20 (right panel). Note that the aspect ratio of the figure is not one. 



changes orientation across the current sheet according to 
By = B tanh(2a;), where the value of Bq is given in terms 
of the magnetization parameter <r m = B^/(2poTq). All over 
the grid there is a small background uniform resistivity r/b, 
except for a circle of radius r v — 0.8, defining a region of 
anomalous resistivity of amplitude rjio = 1.0. As a result, 
the resistivity can be written as 



ri b + T]io [2(r/r v f - 3(r/r^) 2 + l] for r^r v , 



>lb 



for 



r > r„ 



(23) 



where r — *J x 2 + y 2 . The velocity field is initially zero, 
hence To = 1, while the electric field is given by E z = 
Tj{dB y / dx). In most of our simulations we have considered 
the case with po = 1, po = 1. The numerical grid consists 
of an unstructured mesh composed of triangles in 2D and 
tetrahedrons in 3D, which are clustered along the current 
sheet. The grid extension is given by [—50, 50] x [—150, 150] 
in 2D and by [-50,50] x [-75,75] x [-12.5,12.5] in 3D. 
Tab. [T] reports the basic parameters of the models studied 
in our simulations. The Lundquist number S = VAL/rjb is 
reported in the fourth column, where L is the length of the 
initial current sheet, w hile v\ = B 2 /(hp + B 2 ) is the rela- 
tiyistic Alfven velocity l|Komissarovl Il997l ; iDel Zanna et al.1 
2007). The name of the models have been chosen to facili- 
tate recognizing their main parameters. For example, model 
2D-mlO-S2 . 68e5 is a two dimensional model with magneti- 
zation a m = 10 and Lundquist number S = 2.68 x 10 5 (cor- 
responding to r)b — 10 -5 ). We have used periodic boundary 
conditions at y m i n and y max , while zeroth order extrapola- 
tion is applied at :r m i n and x max . 



4.2 Results 

In a first series of simulations we have considered the case 
of an isotropic Ohm law, namely adopting Eq. Q20p with 
f = 0, concentrating on the effects that an increasing mag- 
netization produces on the system. Fig. [3] shows the rest 
mass density and the magnetic field lines for the model 
2D-ml-Sl . 68e5 with a m = 1.0 (left panel) and for the 
model 2D-m20-S2.82e5 with a m = 20 (right panel). The 
two panels, which show snapshots of the two models at 
the same time t = 90, confirm the essential features of a 
Petschek-type relativistic reconnection, as a l so reported by 
other authors (Watanab e fc Yokov ama 2006; Zcn itani et al.l 
l2009bl iah. Namely, magnetic reconnection is triggered in 
the region of anomalous resistivity, producing the typical 
A— type topology of the magnetic field. As a result, mag- 
netic energy is converted into both thermal and kinetic en- 
ergy, and two plasmoids moving in opposite directions, and 
corresponding to the two high-density regions of the figure, 
are accelerated along the direction of the magnetic field. The 
magnetic tension, in fact, is responsible for the collimation of 
the flow. The plasmoid is highly compressed by plasma with 
much higher velocity and lower density (see the discussion 
about Fig. [S] below) and its rest mass density increases with 
increasing magnetization. A complementary information to 
that of Fig. [3] is provided by Fig. [4] The top and the middle 
panels show, respectively, the time evolution of the magnetic 
energy (normalized to its initial value) and of the Lorentz 
factor of the plasmoid, which is also the maximum Lorentz 
factor monitored over the grid. The dissipated magnetic en- 
ergy produces an increase of both the thermal and kinetic 
energy. The latter results in the acceleration of the plasmoid, 
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Figure 4. Time evolution of the magnetic energy (top panel) 
and of the Lorentz factor (middle panel) for the first four mod- 
els reported in Tab. [T] with increasing magnetizations cr m . The 
bottom panel shows the dependence of the Lorentz factor on the 
magnetization in the plane log<r m — logT. 



which is more efficient for higher magnetizations. However, 
as the region around the anomalous resistivity becomes more 
and more rarefied, the conversion of magnetic energy into 
thermal energy becomes more efficient than the conversion 
into kinetic energy and the Lorentz factor reaches a satura- 
tion. The small wiggles visible in the Lorentz factor curves, 
particularly at high magnetizations, are due to impulsive re- 
connection episodes, taking place in very small and localized 
regions around the plasmoid. The bottom panel of Fig. [4] on 
the other hand, shows the maximum Lorentz factor reached 
as a function of the magnetization parameter. Our results 
show that, while at low magnetizations the Lorentz factor 
obeys a dependence of the kind F oc <7 ly/3 , at higher magneti- 
zations the depend ence is close to the theoretical predictio n 
given by F oc all 2 (|Lvubarskvll2005l ; IZenitani et al.ll2009bl h 
The reconnection rate, namely the speed at which the 
reconnection process takes place and that we have computed 
aijf| r = Ez/Bo, is also stro ngly dependent on the magnet i- 
zation, as firstly noticed bv lWatanabe fc Yokova ma (2006). 
We have followed its temporal evolution and found that it 
reaches a maximum around t ~ 20 — 30, while decreasing 
asymptotically after that to reach a stationary value. The 
maximum reconnection rate is r ~ 0.07 and r ~ 0.2 for 



See lZenitani et al ] j2009bh for alternative definitions of the re- 
connection rate. 
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Figure 5. Color map of the rest mass density (left panel) and 
of the Lorentz factor (right panel) for the model 2D-m20-S2 . 82e5 
with a m = 20, S = 2.82 X 10 5 at time t = 90. 

models with a m = 1.0 and a m = 20, respectively. It should 
be noted that the regions of maximum rest mass density 
and of maximum Lorentz factor do not match exactly. The 
mismatch is reported in Fig. [5] which provides a focus of the 
upper plasmoid that is visible in the right panel of Fig. [3] 
The left and the right panels of Fig. [5] show the rest mass 
density and the Lorentz factor, respectively. As it is apparent 
from this figure, the plasma reaches its maximum velocity 
at the basis of the plasmoid in a very rarefied region. On the 
other hand, the portion of the plasmoid with maximum rest 
mass density has very low Lorentz factor. As a result, a tiny 
bow shock is produced at the basis of t he pla smoid, confirm- 
ing similar findings bv lZenitani et all l|2010t ) who performed 
a detailed analysis of the generation of slow shocks around 
the plasmoids. 

In a second series of simulations, we have analyzed the 
dependence of the reconnection process on the background 
Lundquist number 5*, while keeping the same peak value 
of the anomalous resistivity rjio = 1.0 and of the magneti- 
zation parameter <j m = 10. It is worth recalling that high 
values of S, corresponding to higher background conductiv- 
ities, represent a challenge for the numerical scheme, as the 
stiffness of the resistive magnetohydrodynamics equations 
becomes more severe. However, a higher resistivity jump 
between the background and the anomalous resistivity re- 
produces physical conditions where the plasma conductivity 
changes sharply over small distance scales, as expected in re- 
alistic conditions. The top panel of Fig. [6] shows that, when 
increasing the Lundquist number from 5* = 2.68 x 10 5 to 
5* = 2.68 x 10 8 , the asymptotic Lorentz factor reached by 
the plasmoid is smaller. The explanation of this effect is once 
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Figure 6. Time evolution of the Lorentz factor (top panel), of 
the thermal energy (middle panel) and of the reconnection rate 
E z /Bo (bottom panel) for models having the same magnetization 
c m = 10.0 but different Lundquist numbers. 



again to be found in the efficiency of magnetic energy con- 
version. When there is a higher conductivity jump between 
the central hot spot and the background, in fact, the frac- 
tion of magnetic energy that converts into thermal energy 
increases. This is shown in the middle panel of Fig.(S] which 
reports the time evolution of the specific internal energy of 
the plasma, which is higher for higher Lundquist numbers. 
The bottom panel of Fig. [S] on the other hand, shows the 
evolution of the reconnection rate r = E z /Bo, whose peak is 
a mild decreasing function of the Lundquist number. As it is 
evident from the top panel of Fig. [5] the Lorentz factor for 
the model with S = 2.68 x 10 8 manifests an irregular evolu- 
tion after t = 70, which is the signature of the development 
of a magnetic instability. 

This effect is reported in Fig. [7] showing the rest mass 
density in two models having different Lundquist numbers. 
No sign of instability is visible in simulations with Lundquist 
numbers as large as S ~ 10 7 (left panel), while a chain 
of magnetic islands is produced for S > S c ~ 10 8 (right 
panel) . This instability, which resembles a tea ring instability, 
was in vestigated through a linear analysis bv lLoureiro et ail 
(5553), and subsequently confirmed via numerical simu- 
lations in the Newtonian framework by [Sarafancy et al.l 



2009). Our results, combined with those bv lSamtanev et al.l 
2009), indicate that, in the transition to the relativis- 
tic regime, the critical Lundquist number increases from 
S c > 10 4 to S c > 10 8 . Such a conclusion may have deep 
implications for high Lundquist number reconnection in re- 
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Figure 7. Generation of plasmoid chains in very high Lundquist 
numbers configurations. Left and right panels report, respectively, 
the color map of the rest mass density at time t = 60 for two 
models having S = 2.68 X 10 7 and S = 2.68 X 10 s . 



alistic astrophysical conditions, and it will deserve further 
investigations. 

Finally, in a third series of simulations we have consid- 
ered the effects of the full anisotropic Ohm law given by 
Eq. (|20p with £ 7^ 0. In this case, a component of the mag- 
netic field perpendicular to the x — y plane must be intro- 
duced, otherwise the term E-B in (|20[) remains zero. Hence, 
we have investigated anisotropic effects in three spatial di- 
mensions and in configuratio ns with B z == , the so called 
"guide field" configurations (|Lvubarskv||2005D . for which we 
have chosen the value B z = 0.5Bo. As a representative ex- 
ample, Fig.[8]shows the magnetic field B 2 on the slices x = 
and z — for the model 3D-ml . 25-BzO . 5 at time t — 90. 
The diffusion of the magnetic field is clearly visible in the re- 
gion around the anomalous resistivity. It has to be remarked 
that, within our single fluid approxi mation, inertial r esistiv - 
ity effects like those encountered bv lZenitani et all (|2009a[ ) 
cannot be captured and are not discussed here. Simulations 
with the guide field configuration in combination with the 
anisotropic Ohm law turned out to be very challenging for 
the numerical scheme. In particular, strong magnetizations 
could not be reached and the results that we show here are 
limited to the case o m ~ 1. When discussing the effects in- 
troduced by an anisotropic Ohm law, we first need to distin- 
guish them from thos e produced by the guid e-field. By using 
a two-fluid approach. IZenitani et al.l l|200 9a) concluded that 
the guide field modifies the composition of the output en- 
ergy flux, which, as the guide field increases, changes from 
being enthalpy dominated to be Poynting dominated. Even 
within a single fluid approach, we confirm here a similar 
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Figure 8. Magnetic field B 2 on the slice x = and on the slice 
z = for the model 3D-ml . 25-BzO . 5 at time t = 90. 



finding. This is shown in Fig. [9] where we compare the be- 
havior of three different configurations each of which with 
an effective cr m = 1.25. The solid red and the dashed blue 
line, in fact, refer, respectively, to a 3D simulation without 
the guide magnetic field and to a 3D simulation with the 
guide magnetic field, but with an otherwise isotropic Ohm 
law. In the case of a non zero guide field, the Lorentz fac- 
tor is substantially smaller, while the decay of the magnetic 
energy, not reported here, is corre spondingly slower. Th is 
is consistent with what reported bv lZenitani et all (l2009al ). 
who found that the introduction of the guide field has the 
net effect of diminishing the bilk kinetic energy. Having clar- 
ified this, we have increased the parameter £ to establish the 
extent to which results are affected by anisotropic Ohm law 
(green long-dashed line) . The present implementation of the 
numerical scheme does not allow to treat the anisotropic 
parameter £ as a true free parameter, and we have there- 
fore concentrated on a mildly anisotropic regime. When £ 2 
is increased from (isotropic case) to 0.5, the Lorentz factor 
indeed grows, but by 1% only after t = 100. Intuitively, this 
effect can be explained in terms of the increased Lorentz 
force on the plasma, because of the extra current term on 
the right hand side of (|2ip . More work is needed to analyze 
the anisotropic regime under more realistic conditions. 



5 CONCLUSIONS 

We have investigated the dynamics of High Lundquist num- 
ber relativistic magnetic reconnection, by performing two 
dimensional and three dimensional magnetohydrodynamics 
simulations of a Petschek type configuration. By resorting 
to high order discontinuous G alerkin methods as proposed 
bv lDumbser fc Zanottil (2009), we have found that Lorentz 
factors up to ~ 4 can be obtained for plasma parameters 
<7 m up to 20 in systems with Lundquist number S ~ 10 5 . 
However, such values of the Lorentz factor decrease when 
S is increased as a result of a reduced background resistiv- 
ity. When S is larger than a critical value S c ~ 10 8 , the 
Sweet-Parker layer becomes unstable, generating a chain of 



Figure 9. Combined effects of the guide field and of anisotropic 
Ohm law on the time evolution of the Lorentz factor. All of the 
models shown have a m = 1.25. See text for explanation. 



secondary magnetic islands. This effect deserves additional 
theoretical investigations, in view of the fact that a simi- 
lar instabili ty has already been r eported in the Newtonian 
framework l|Samtanev et al.| [2009). but at much lower criti- 
cal Lundquist numbers, i.e. for S > S c ~ 10 4 . 

We have also shown that, when an anisotropic Ohm 
law is adopted, both the reconnection rate and the Lorentz 
factor of the accelerated plasmoid are slightly increased. 
Although this increase is small and within a few percent 
with respect to the isotropic Ohm law, strongly anisotropic 
regimes remain challenging even for our advanced numer- 
ical scheme. We plan to improve on these limitations in 
our future applications of relativistic magnetic reconnection 
to the curved spacetime around a neutron star, where an 
anisotropic Ohm law can play a substantial role. 

Finally, it is worth mentioning that pure magne- 
tohydrodynamics reconnection occurs in collisional plas- 
mas, while the transition to collisionless re connection 
is likely to enhance the reconnection rate (|Uzdenskvl 
120071; lUzdenskv et"ai1 |201Ch . and, presumably, also the 
acceleration properties. This possibility, whic h is related 

to so me recent results of plasma physics l|Hesse et alj 

20091; iBessho fe Bhattacharieel 120071 ) . has been explored by 



Goodman fc Uzdenskvl ll200ct) in a ccretion disc coronae and 
bv iMcKinnev fc Uzdenskvl d2010h to trigger dissipation in 
Gamma-Ray Burst jets. The numerical exploration of colli- 
sionless reconnection has received less attention^ and it will 
represent another direction of our future analysis. 



3 Also note that Lazarian et al. (2010), after solving numerically 
a reduced set of equations, showed that reconnection in a turbu- 
lent fluid occurs at a rate comparable to the rms velocity of the 
turbulence, irrespective of the degree of collisionality. 
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APPENDIX A: DERIVATION OF EQ. (21) 

We first write the fluid four velocity, the magnetic field in 
the comoving frame and the electric field in the comoving 
frame as 

u" = r(/+/) (Al) 
gM = p ^ Uv ^Y^^.^ + E^ + iHx BY] (A2) 
fe" = F^Uu = T [n M (v ■ B) + 5" - (v x EY] (A3) 

where we have used the definitions $Jfy and ((3} for the elec- 
tromagnetic tensor and where we have used the property 
<^ vX = e M " AK n«. We then refer to the relation 

a ^ = a ( g ^ + fbW) (A4) 

and after noticing that 6 M e M = E ■ B (it is a relativistic 
invariant) we rewrite (|18[l as 

F = q T(v t ' + n M ) + aT[(v ■ £)n M + £ M + (v x B) M ] + 
af(v ■ E)T[(v- B)n' 1 + - (v x £)"] . (A5) 

At this point we compare (|A5[) with J M = p c n M + J* 1 to find 

Pc = qoT + aF(v ■ E) + a^ 2 F(E ■ B)(v ■ B) (A6) 
J M = q Tv» + oTE^ + aF(v X B) M + 

<T£ 2 {E-B)r(B> 1 -{vxEY) ■ (A7) 
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A last replacement of go from (|A6|) into (|A7|) allows to derive 
(|20|) of the main text. 



